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We study the dynamics of small fluctuations about the uni- 
form state of a crystal moving through a dissipative medium, 
e.g. a sedimenting colloidal crystal or a moving flux lattice, 
using a set of continuum equations for the displacement fields, 
and a one-dimensional driven lattice-gas model for the cou- 
pled concentration and tilt fields. For the colloidal crystal 
we predict a continuous nonequilibrium phase transition to a 
clumped state above a critical Peclet number. 
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What is the response of a crystalline lattice to a small, 
long-wavelength disturbance? Elastic theory jjj and its 
extension to time-dependent phenomena [|2j provide a 
complete answer to this question for a system in ther- 
mal equilibrium. In this Letter we ask and answer the 
same question for a lattice being driven through a dis- 
sipative medium by a constant external force. This im- 
portant nonequilibrium steady state arises, for example, 
in the steady sedimentation ||[| of a colloidal crystal || 
and in the motion of a depinned flux lattice in a current- 
carrying superconductor. Using symmetry arguments we 
construct continuum and lattice-gas models for the dy- 
namics of small distortions about the uniform state of 
such a steadily moving lattice. The effects that we dis- 
cuss arise from the dependence of the mobility of a given 
region on the local strain of the crystal. Our most strik- 
ing result is that the dominant linear response at long 
wavelengths is proportional to the driving speed of the 
lattice, not to its elastic constants, and that this response 
can lead to a nonequilibrium phase transition. 

Before presenting our results in more detail let us re- 
call an important early study. Crowley || carried out 
experiments (on steel balls dropped gently into turpen- 
tine oil) and theory (calculating the hydrodynamic in- 
teractions between the spheres) to show that a regular 
horizontal array of sedimenting spheres was linearly un- 
stable to clumping and buckling. Elastic forces, Brownian 
motion, and nonlinearities, all of which can resist this in- 
stability, were not considered in || . Since experiments on 
model colloidal systems are most conveniently performed 
by varying interaction strength rather than temperature 
|^| , the appropriate Peclet number for this problem is 
the ratio Pe — Apag/G of gravitational to elastic rather 
than Brownian stresses. Here Ap is the difference be- 
tween the mass-densities of particles and solvent, a the 
particle radius, g the acceleration due to gravity, and G a 
typical elastic modulus of the suspension. A sedimenting 



colloidal crystal, according to || , is therefore unstable in 
the Pe = oo limit. Hence the question posed in our title: 
does the instability persist at finite Pe? 

Our answer to this question is in two parts: (i) analyt- 
ical results from a system of coupled nonlinear stochas- 
tic partial differential equations for displacement fields, 
and (ii) numerical studies of an equivalent driven lattice 
gas model M describing the coupled dynamics of con- 
centration and tilt fields. Our simulations of the lattice 
model are done in the parameter range where there is 
a local tendency to undergo the Crowley || instability. 
We find that the system is nonetheless stable with re- 
spect to clumping upto a critical Pe at which a contin- 
uous nonequilibrium phase transition to a clumped state 
takes place. For the corresponding parameter range the 
continuum model, neglecting nonlinearities and noise, is 
unstable without threshold, i.e. at Pe = 0. In view of 
the results of the lattice model, we expect that nonlinear- 
ities and fluctuations shift the critical Pe for the onset of 
clumping in the continuum model to a nonzero value. 

We first derive the continuum model and perform a 
mean- field analysis, then show how the lattice-gas model 
follows naturally from it. Ignoring inertia altogether, 
which is justified for the experimental geometry we wish 
to consider |11|, the displacement field u(r, t) of a d- 
dimensional lattice moving through a frictional medium 
with a mobility which depends on the local strain obeys 



ii = /x(Vu)(DVVu + F + C) 



(1) 



where the first term on the right represents, through the 
tensor D, the elastic restoring forces, the second is the 
driving force, and the third is a random force of thermal 
or possibly hydrodynamic origin [p"2f . fi is the mobility 
tensor which, in the single particle limit for colloids in a 
solvent with viscosity n, approaches the isotropic Stokes's 
Law value 6 * ■ We expand fi in powers of Vu: 



M(Vu) = Mo + AVu + C((Vu) 2 ), 



(2) 



where fio is the mean macroscopic mobility of the undis- 
torted crystal. For steady sedimentation along z, assum- 
ing isotropy in the d — 1 transverse (_L) dimensions but 
not under z — > — z, (|l|) and (||) lead directly to 

u'j_ = \id z u± + A 2 V_lu z 

+0(VVw) + 0(VuVu) + fL, (3a) 

Uz = A3V±.U_l + \id z U z 

+C(VVu) + O(VuVu) + f z , (3b) 
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where the constant drift along z has been removed by co- 
moving with the crystal |1 . Here f is a spatiotemporally 
white noise source |tJ and \ are phenomenological coeffi- 
cients whose origin, in the case of a colloidal crystal, lies 
in the hydrodynamic interaction between the particles 

filp^] . We explain the physical content of the terms in 
below. First note that linearizing and Fourier trans- 
forming in space and time yields modes with frequencies 
of the form uj = ±\ / X 2 Xsq — iDq 2 for q z = 0, q± — > 0. 
If A2A3 > 0, this leads to wavelike excitations at small 
wavenumber q, which are not the usual shear-waves of 
a crystal at equilibrium. The latter have been turned 
already into diffusive modes by the frictional dynamics 
adopted in (||). The speed of these waves is determined 
by the strain-dependence of the mobility, and the damp- 
ing by the tensor D which is a ratio of elastic constants 
to friction coefficients. 

When A2A3 < 0, the dispersion relation at small q 
becomes oj ~ ±iq so that the model is linearly unstable 
with growth rate oc q. Thus there exist long-wavelength 
distortions of the perfect lattice which grow exponentially 
in time within the linear theory. We do not know the sign 
of A2A3 for moving flux lattices, but for colloidal crystals 
it is negative, making them linearly unstable. This is 
because hydrodynamic interactions cause denser regions 
in the suspension to sink faster, and tilted regions to 
move laterally in a way so as to cause an instability ||. 
Note that in ^ the linear elasticity of the crystal enters 
at second order in wavenumber and can thus not alter our 
conclusions about linear stability at long wavelengths. 
For q > ~ VA2A3/D, elastic forces suppress the linear 
instability. Small crystals are thus linearly stable. 

To go beyond this linear analysis is daunting: even in 
d = 2, symmetry permits nine terms bilinear in Vu and 
six linear second derivative terms. Remarkably, all the 
essential physics is retained in a greatly simplified ver- 
sion in one space dimension. Consider a two-dimensional 
crystal described by a two-component displacement field 
{u x , u z ), with the sedimentation direction z averaged out 
so that only x derivatives are included. The equations 
of motion, retaining the lowest order nonlinearities and 
gradients, read 



Ux = X 2 d x u z + jid x u x d x u z 
u z = X 3 d x u x + j 2 (d x u x ) 2 
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+l3,{d x u z 



D 2 d x 



Did x 2 u x + f x 



(4a) 
(4b) 



The physics of each of the terms in (||) or (Q) is rea- 
sonably clear. The first two terms on the right of ( f4a| ) 
say that a tilt (i.e. d x u z ) produces a lateral drift. The 
first two terms on the right of (^) contain the concen- 
tration dependence, and the third (a Burgers/KPZ-like 
term |l(J]) the tilt-dependence, of the vertical sedimenta- 
tion speed. Note that the 7, terms can be seen as arising 
from the dependence of the A^s on concentration and tilt 
. The second derivative terms in both equations arise 



simply from elastic forces, and f x and f z are spatiotem- 
porally white noises. The symmetry of (Q) is under the 
joint inversion x — > — x, u x — > —u x . 

For ~ jj, u x — ► u x — X 2 xj^\ eliminates the 
A2 and A3 terms from (||), reducing it to the model of 
Erta§ and Kardar (EK) (l4) in their limit Aj^ = An, 
with its higher symmetry (x — > —x), albeit in a state 
of nonzero mean d x u x . The fluctuation-dissipation the- 
orem, Galilean invariance and Cole-Hopf properties that 
arise in Jl4| for special parameter values thus obtain 
here as well. If in addition 71 = 27 3 and D\ = D 2 , 
the equations decouple in terms of transformed vari- 



ables 



u x ± 



into two equations: 



±V\2\3d x (/)±+D 1 d x <j>±±yj ^72(9 X (}>±) +f±, a pair of 
independent KPZ [1(| equations with oppositely directed 
kinematic wave |lj ] terms, and nonlinear couplings with 
opposite signs. Clearly, for these parameter values, the 
stable driven crystal should exhibit KPZ exponents in 
its correlation functions. The relevance of perturbations 
about the highly symmetric EK limit, as well as the stat- 
ics and dynamics of the 'stable' case A2A3 > will be 
studied in later work. In the present paper we focus on 
A 2 A 3 < 0. 

Let us first look for steady-state solutions to (Q) in 
the absence of noise, in terms of p = d x u x (the local 
concentration fluctuation) and 9 = d x u z (the local up 
or down tilt). If we restrict ourselves for simplicity to 
spatially uniform states with left-right symmetry (so that 
9 = 0, and the net currents of p and 9 are zero), we 
are left with only two possibilities: p = 9 = 0ovp = 
— A3/72, 9 = 0. In the vicinity of r = A2A3 = 0, the 
former is stable for r > 0, the latter for r < 0. This 
exchange of stabilities leads to a continuous onset of the 
p 7^ state, p ~ |rp with [3=1. Similar analysis |l6| ] 
gives a correlation length diverging as \r\~ v with v = 1/2. 

Instead of attempting a perturbative treatment of the 
effect of nonlinearities and fluctuations on the above 
mean-field picture, we replace the continuous variable x 
by a discrete index i, and du x /dx by pi = u x (i+l)—u x (i) 
[similarly 9i = u z (i + l) — u z (i)], with pi and 9i restricted 
to ±1. Such an approach || has proved very success- 
ful for simulating the KPZ equation. The "paramag- 
netic" phase of these Ising variables corresponds to the 
the undistorted crystal, and the "ferromagnetic" phase 
represents a macroscopically clumped and tilted state, in 
terms of suitable order parameters which we define below. 
The best way to visualize the discrete model is to think of 
two sublattices: a typical configuration can then be de- 
scribed by a sequence of spins p\9\p 2 92Pz9z . . .. The dy- 
namics of the spins is constructed by analogy with lattice 
models [|| for the KPZ equation, retaining the essential 
features of (Q), viz., conservation of 9 and p, stochasticity, 
lack of up-down symmetry, and the bias provided by each 
species on the motion of the other. The two approaches 
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should yield identical long-distance properties. 

Let us denote the states of pi by '+','— ' and those of Bi 
by '/' (uptilt) and '\' (downtilt). In the update rule cor- 
responding to the linearly unstable case of (|]), the rates 
for the following exchanges are enhanced relative to the 

corresponding reverse rates: +\ ► — \+, — /-\ > +/—, 

I + \ — ► \ + / and \ — /—*/ — \ Since we are modelling 
charge-stabilised suspensions, it is useful to introduce a 
repulsion between regions of high density in the form of 

an enhanced probability for + H — or h + to go to 

H K Combining all of the above, we get the following 

exchange probabilities for adjacent pairs of concentration 
and tilt: 



P 



Pi,Pi+i 



D n 



9iPi + a{(l + pi)(l + 



Po. 



+ (l + p i+ i)(l+p i+2 )}, 
4+ i — D g + e s pi+i9i + g\Pi + \ + gzQi, 



(5a) 
(5b) 



where D p and Dg are related to the elastic constants, e p , 
eg, gi and g2 to the 7^ and Xi in (^), and a is the repulsion 
[|l7[| . Note that decreasing a reduces the stiffness of the 
system, thus increasing the effective Peclet number. The 
last two terms in (^) arise because of the lack of up-down 
symmetry. Our results in this paper are for e p eg > 0, 
corresponding to A2A3 < in (Q). 

The mean value of both p and 6 are expected to be 
zero in the experimental system. We worked, therefore, 
at zero total 'magnetization' for both fields and stud- 
ied the model starting from random initial conditions, 
evolving it according to the above update rules for vari- 
ous system sizes TV. Periodic boundary conditions were 
used for all the runs. For high values of repulsion the 
spin configurations continued to be homogeneous under 
time evolution. When the repulsion was small or absent, 
there was a phase separation into regions of high and low 
concentration and of up and down tilt, separated by in- 
terfaces. Thus the lattice seems to be stable for strong 
repulsion, but undergoes Crowley's clumping instability 
|^| for weak repulsion. The same behaviour, qualitatively, 
is observed as D or Dg are increased keeping other pa- 
rameters fixed fllql . 

To describe the ordered phase in this model with con- 
served dynamics, we use two essentially equivalent or- 
der parameters: $ p = J ~jfJ2iPiPi+N/2 (similarly $ e 
for tilt) which measure how anticorrelated the spins are 
across half the system size; and the moduli 

of the Fourier amplitudes of the spin fields at the smallest 
nonzero wavevector k± = ^ jl8[ (the amplitude at k = 
is zero). Fig. |l| shows that the order parameter <J> p is 
appreciable for small repulsion and decreases rapidly to 
a value consistent with zero for sufficiently large repul- 
sion a. Moreover, <& p increases with N for a small and 



decreases 16 roughly as X/\N for a large. There must 
thus be a continuous nonequilibrium phase transition at 
a around .05, although to pin down the critical value of 
a would require careful finite-size scaling. 
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FIG. 1. Order parameter as a function of repulsion 

strength a for various system sizes N. Increasing a is like 
decreasing Peclet number. 



We now present an independent check that the ob- 
served phase-separation is not merely the result of tran- 
sients. A truly phase-separated state in a system of 
length N should have barriers to remixing which grow 
as for some power £. The lifetime of such a state 
would then go as exp(iV>). To look for such barriers, we 
define a lifetime t(N) to be the mean time of first pas- 
sage of the order parameter from a value 6 2 |*| max 
to &i l^lmax, where l^lmax is the maximum value of 1^1 
over the runs. b\ and 62 are numbers independent of N, 
chosen to get good statistics; we took b 2 — 0.8, b x = 0.5. 

Fig. H shows that t(N) is roughly exponential in N for 
a = and distinctly faster than a power law, presumably 
a stretched exponential, for a = 0.015. This is strong 
evidence ]l9[ for a true clumped phase at weak repulsion. 

To see why phase separation can occur in this one- 
dimensional model one has to look at the positions of 
the concentration and tilt domains. We find in our simu- 
lations that the system goes into a steady state in which 
the domains are staggered with respect to each other by 
an approximate distance N/A. This happens in such a 

way that a concentration interface + + H lives in a 

region crowded with uptilts / which inhibit the exchange 
of a pair H — . The dissolution of the interface by inter- 
diffusion of + and — thus requires uphill motion over a 
nonzero fraction of N. 

Since our simple one-dimensional model undergoes a 
clumping transition, it is reasonable to expect that a 
real charge-stabilized colloidal crystal in a fluidized bed 
0,0 will do so as well. The repulsion between poly- 
balls may be decreased by adding salt to the fluid, which 
should lead to an observable clumping transition at ionic 
strengths much lower than those required to produce 
melting or aggregation at equilibrium. The clumping will 
manifest itself as a breakup of the crystal into smaller 
crystallites (since the crystal is stable at small enough 
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Lifetime of the ordered phase r(N), from the 
concentration order parameter \^ p \. a = 0: (a) log-log and 
(b) semi-log plot, a = 0.015: (c) log-log and (d) semi-log 
plot. Similar results hold for the tilt order parameter Jl6{ . 
Note the upward curvature of the log-log plots in either case, 
indicating lifetimes growing faster than any power of TV. The 
downward curvature in (d) indicates that r(N) is slower than 
for a — .015. The error in r is determined by allowing the 
system to run for several lifetimes and statistically estimating 
the standard deviation. 



system size) , separated by regions of strong upward fluid 
now f20[| . A detailed analysis of this dynamics would 
require the inclusion of the hydrodynamic flow. 

In summary, we have demonstrated that the long wave- 
length dynamics of a crystal moving steadily through a 
dissipative medium is qualitatively different from its equi- 
librium counterpart. In particular, we have shown that a 
natural driven lattice gas model for this system shows a 
dramatic nonequilibrium phase transition to a clumped 
state, and we urge experimenters to test our predictions. 
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